'''
Author: Zhuangyu
Date: 2022-06-13 16:50:36
LastEditTime: 2022-06-13 16:50:36
'''
'''
Author: Zhuangyu
Date: 2022-06-02 14:31:11
LastEditTime: 2022-06-13 16:50:03
'''
from __future__ import print_function, division
import time
#Importing relevant modules and packages
import numpy as np
import matplotlib.pyplot as plt
import random
from RichardsModel_1D import *
from Parameters import *
from scipy import integrate, linalg, interpolate
from tqdm import tqdm
import mpctools as mpc
import casadi
import scipy.io as sio  

p=Loam() #Soil parameters

#Space and time parameters
Nr,Nt,Nz,dr,dt,dz,Np,NN,Nx,Nw,Ny,Nv,Hz=circular_parameters()
DeltaT,Nsim=time_parameters()

Kc=0.88*np.ones(Nsim)      #Crop coefficient

Et=1.389e-08*np.ones(Nsim) #Reference evapotranspiration

def ode(x,t,u,w,ui,Kc,ET0,p):  
    return RichardsModel_1D(x,u,w,ui,Kc,ET0,p)



u=np.zeros(6*120)
uu=np.zeros(Nsim)
for i in range(6*120):
    if i < 8:
        u[i]=-7e-06 #Moisture is supplied only at the first time instant
    else:
        u[i]=0
uu=np.tile(u,(1,60))
uu=np.transpose(uu)


#Initial condition [pressure head]
x0=-0.8*np.ones(Nz)

#Numpy array to store the simulation results [pressure head]
xx=np.zeros((Nsim+1,Nz))

ode1=np.zeros((Nsim+1,Nz))

#Model uncertainty

np.random.seed(1)
ww=0.0000*np.random.randn(Nsim,Nz) #If no model uncertainty is present, ww is multiplied by 0
np.random.seed(2)
www=0.000001*np.random.randn(Nsim,Nz)
np.random.seed(3)
vv=0.01*np.random.randn(Nsim+1,Nz) #measurement noise

#Numpy array to store the simulation results [volumetric moisture content]
yy=np.zeros((Nsim+1,Nz))
yy[0,:]=getOutputs1D_allnodes(x0,p).ravel()+vv[0,:]

#unknown inputs
a1=0*np.ones((Nsim+1,Nz))
a=0.00003*np.ones((Nsim+1,Nz))

e0 = time.time()

for i in tqdm(range(1, Nsim+1)):
    # Generate observed data (x bias)
    ode1[i-1,:]= RichardsModel_1D(tuple(xx[i-1]), uu[i-1],www[i-1],a1[i-1],Kc[i-1],Et[i-1],p)
    xx[i,:]=xx[i-1,:]+DeltaT*ode1[i-1,:]+a[i-1,:]
    yy[i,:]=getOutputs1D_allnodes(xx[i,:],p).ravel()+vv[i,:]
    
t= list(range(Nsim+1))
# plt.figure(1)
# for i in range (Nz):
#     plt.subplot(5, 7, i+1)
#     plt.plot(t,xx[:,i],'r',t,xx1[:,i],'b')
    

plt.figure(2)
for i in range (Nz):
    plt.subplot(5, 7, i+1)
    plt.plot(t,xx[:,i],'r')
plt.show()


